Series de Tiempo I–Notas de Curso

Universidad Autónoma Metropolitana

email: mperez.mendez@izt.uam.mx

url: https://github.com/MarcoAPerezM

url_lab: https://sites.google.com/view/deserlab/

1 Introducción

El análisis de series de tiempo se configura como una herramienta fundamental entre el conjunto de herramientas con las que cuenta el economista moderno. El mercado laboral del economista demanda el uso de instrumentos analíticos que le permita realizar un diagnóstico preciso sobre los diferentes fenómenos económicos, políticos, sociales y financieros. En los tiempos modernos, las herramientas se enfocan, cada vez más, en un uso intensivo de datos acompañado de capacidades analíticas derivadas de la aplicación de la teoría económica. Entre estas se encuentra el modelado matemático de fenómenos y su correspondiente contrastación empírica. El análisis de series de tiempo está, entre las herramientas econométricas, ubicado como un conjunto de instrumentos que el mercado laboral demanda y valora en demasía. En este libro se encuentra el desarrollo de las series de tiempo univariadas. En el capítulo 1 se desarrollan las características de las series de tiempo, las herramientas básicas para su análisis y se esgrime la necesidad de desarrollar modelos sofisticados, pero sencillos, que permitan entender el pasado de los fenómenos y predecir su futuro. En el capítulo 2 se desarrolla el Modelo de Componentes Subyacentes: tendencia, ciclo y estacionalidad se modelan desde la perspectiva clásica. En el capítulo 3 se desarrollan filtros y suavizamientos. Medias móviles simples, exponenciales, etc. muestran como se comportan los fenómenos en el largo plazo. Se muestra como el filtro Holt-Winters permite reproducir con un alto grado de precisión los fenómenos con tendencia y estacionalidad y, además, permite obtener pronósticos altamente eficientes. En el capítulo 4 se desarrolla la metodología Box-Jenkins y se construyen los modelos tipo SARIMA.

1.1 Requisitos

Como es bien sabido, todo Análisis ecónomico de series de tiempo requiere de ciertos prerequisitos imprescindibles para poder desarollar el análisis en cuestión. En particular, este libro requiere el conocimiento y dominio de diferentes tópicos que le son imprescindibles al economista moderno.

En primera instancia, se encuentra el entendimiento de herramientas estadísticas y de probabilidad que se requieren, tanto para el desarrollo de instrumentos como de indicadores que permitan reproducir los patrones de comportamiento subyacente de los fenómenos económicos, al mismo tiempo sirven para evaluar el grado de exactitud de sus reproducciones.

También se requiere el entedimiento, en un nivel básico, de teoría econométrica: estimación puntual, por intervalos, por cuantiles, etc. Técnicas de estimación de parámetros y modelado matemático.

1.2 Datos económicos y financieros

Otro de los requisitos fundamentales es saber identificar con precisión las fuentes de datos donde se encuentran las series de tiempo relevantes para el análisis económico.

Cuando el profesional cuenta con datos, tiene dos opciones: 1) la primera es cargar directamente en formato csv (Variables Separadas por Comas) y la segunda es cargarglos en formato xlsx (Excel).

Para abrir datos en csv se debe emplear la función read.csv:

#datos <- read.csv("dirección del archivo/archivo.csv, sep=";")
#datos

Para abrir datos en Excel se debe emplear la función read_excel() perteneciente a la librería “readxl”, así como indicar la hoja en la que se encuentran los datos:

# install.packages("readxl")
#library(readxl)
#datos2 <-read_excel("/dirección/libro.xlsx", sheet="hoja")
#datos2

Una manera mucho mas eficiente y moderna de obtener datos es descargarlos directamente desde repositorios oficiales en donde se encuentran los datos. Por ejemplo, la libreria Quantmod (https://www.quantmod.com. Quantitative Financial Modelling and Trading Framework for R) brinda acceso a las bases de datos de Google, Yahoo, oanda y la FRED (Federal Reserve of Economic Data).

Por ejemplo, por medio de la librería getFX() es posible descargar el tipo de cambio del peso mexicano contra el dolar estadounidense desde oanda:

library(quantmod)

## Loading required package: xts

## Warning: package 'xts' was built under R version 4.3.3

## Loading required package: zoo

## 
## Attaching package: 'zoo'

## The following objects are masked from 'package:base':
## 
##     as.Date, as.Date.numeric

## Loading required package: TTR

## Registered S3 method overwritten by 'quantmod':
##   method            from
##   as.zoo.data.frame zoo

getFX("USD/MXN")

## [1] "USD/MXN"

plot(USDMXN)

o el rublo ruso contra el dolar estadounidense:

getFX("USD/RUB")

## [1] "USD/RUB"

plot(USDRUB)

o el precio del oro en relación con el dolar estadounidense:

getFX("XAU/USD")

## [1] "XAU/USD"

plot(XAUUSD)

Con la función getSymbols() se puede desarcargar el Bitcoin desde la FRED, sólo basta con identificar la clave de la variable de interés, en este caso “CBBTCUSD”:

getSymbols("CBBTCUSD", src = "FRED") #source = src

## [1] "CBBTCUSD"

plot(CBBTCUSD)

1.3 Descargar datos de INEGI

En este caso, el interés principal se ecuentra en analizar series de tiempo mexicanas. Para ello, es necesario conectarse a la API de INEGI (Instituto Nacional de Geografía e Informática https://www.inegi.org.mx/servicios/api_indicadores.html). Una API (Interfaz de Programación de Aplicaciones) es, en términos generales, un conjunto de protocolos o reglas que permiten la interacción entre diferentes aplicaciones de software con el objetivo de que se comuniquen entre sí para intercambiar, principalmente, datos, características y funcionalidades. Para tener acceso a la API de INEGI es necesario solicitar una llave de acceso conocida como token.

Existe una librería que permite la conexión entre R y la API de INEGI.

library(inegiR)
#token <- "token de INEGI"

Si es de nuestro interés analizar el comportamiento del Poducto Interno Bruto Trimestral, es necesario identificar su ubicación en el Banco de Información Económica. La ubicación exacta es Indicadores económicos de coyuntura > Producto interno bruto trimestral, base 2018 > Series Originales > Valores a precios de 2018. EL identificador del PIB trimestral es “735879”. La función para realizar esta conexión es inegi_series(). Es necesario indicar el identificador de la variable, el token y la base de datos, la cual puede ser el Banco de Indicadores o el Banco de Información Económica.

inegi_id <- "735879"
pib <- inegi_series(inegi_id, token, database = "BIE")

Como puede observarse, por medio de la función plot(), la descarga se realiza con el orden inverso, la tabla se encuentra ordenada de los valores más actuales a los más antiguos.

plot(pib$values, type="l")

Es necesario verificar la primer fecha de publicación de la variable y la frecuencia con la que se repite dentro de un año.

tail(pib)

##           date date_shortcut   values notes
## 174 1981-04-01            Q2 11564025    NA
## 175 1981-01-01            Q1 11345848    NA
## 176 1980-10-01            Q4 10927666    NA
## 177 1980-07-01            Q3 10392733    NA
## 178 1980-04-01            Q2 10342350    NA
## 179 1980-01-01            Q1 10401368    NA

Como puede observarse, la fecha de inicio es el primer trimestre de 1980. Es necesario invertir el orden, para ello diseñamos una función que cambie el orden de acuerdo a la fecha y mantenga solo la columna de valores. A esta función personalizada la hemos llamado ts_invert().

ts_invert <- function(x){
  x <- x[order(x$date),]
  x <- x$values
}

Al aplicar esta función personalzida sobre la variable pib tenemos el orden adecuado

pib <- ts_invert(pib)
plot(pib, type="l")

Ahora es necesario transformar la variable en una serie de tiempo. Para ello usamos la función ts() la cual permite identificar la fecha de inico y la frecuencia. En este caso, la variable pib se transformará en una serie de tiempo que inicia en el primer trimestre de 1980 y tiene frecuencia 4, ya que es trimestral.

pib <- ts(pib, frequency = 4, start =c(1980,1) )
ts.plot(pib)

1.4 Indicadores de precisión

Al momento de realizasr un análisis de series de tiempo siempre se va a cometer algún tipo de error. El más frecuente es el error de predicción o de estimación, este se define como la diferencia entre el valor observado y el valor estimado:

\[\hat{e}_t = y_{t} − \hat{y}_t \]

Estos errores son causados por un gran número de fuentes, por ejemplo, hay errores de medición en los procesos de recolección de las variables, hay errores de muestreo, hay errores derivados de la aleatoriedad del comportamiento humano, etc. Estos errores pueden ser grandes y sistemáticos y esto depende, en gran medida de:

a) Identificación errónea de los patrones subyacentes y relaciones entre variables. Podemos identificar, de manera equivocada, un patrón equivodado. También podemos construir un modelo estadístico basado en pocos datos y suponer que el patrón es estable a lo largo de mucho tiempo.

b) Patrones inexactos. En general, en las ciencias sociales y, en particular, en la Economía los patronos que podemos identificar son, por naturaleza, inexactos e imprecisos. Esto se debe a que no es una ciencia determinista como la física o la quimica. Las relaciones económicas no son deterministas, son aleatorias y a pesar de que podemos identificar patrones de comportamiento estables en el tiempo, siempre habrá incertidumbre y múltiples factores que generan variaciones en los resultados. Por ejemplo, podemos construir un modelo que capture, relativamente bien, el ingreso de los hogares en México. Sin embargo, ¿cuántas variables se encuentran involucradas en el la determinación del ingreso de un hogar particular? Supongamos que el ingreso promedio de un hogar es de $12,000. ¿Cuántos hogares ganan esa cantidad? De dichos hogares, ¿cuántos tiene jefatura femenina? ¿Cuántos jefes de hogar tienen educación universitaria? etc. etc. Hay un gran número de variables involucradas en la determinación del ingreso, todas esas variables generan una fuerte variabilidad e impiden que el fenómeno sea determinista, por ello es imposible identificar con exactitud el ingreso de los hogares con base en un modelo estadístico o econométrico. Esta variación induce un error de estimación.

c) Cambios en los patrones. Siguiendo el ejemplo del ingreso de los hogares, un error común es suponer que los patrones que se identificaron en un momento en el tiempo serán los mismos durante un periodo de largo plazo. Esto no es necesariamente cierto, hay cambios en la estructura de los hogares, en la estructura productivia, en los niveles de educación y productividad, etc. Estos cambios inducen cambios en los patrones subyacentes de los fenómenos y por lo tanto en los parámetros de los modelos, a este dilema normalmente se le conoce como cambio estructural o paramétrico.

Por ello es necesario desarrollar herramientas que nos permitan identificar el grado de precisión de los modelos construidos. Para evaluar la precisión vamos a suponer que tenemos un fenómeno \(y_t\) observado.

y <- c(25,28,27,30,35,44,37,41,45)
plot(y, type="l")

Vamos a suponer que tenemos dos modelos alternativos para reproducir y, eventualmente, predecir el fenómeno \(y_t\).

yp1 <- c(23,25,26,29,31,34,41,40,38,44)
yp2 <- c(24,26,28,29,33,34,45,38,43,45)

y <- ts(y)
yp1 <- ts(yp1)
yp2 <- ts(yp2)

ts.plot(y,yp1,yp2, col=1:3)

error1 <- y-yp1
error2 <- y-yp2

1.4.1 Error Medio Absoluto

\[ EMA =\frac{1}{T} \Sigma_{t=1}^T | e_t |\] El EMA es una medida de precisión general que brinda una idea del grado de dispersión y cuenta con la característica de que le brinda la misma importancia a todos los errores.

ema_1 <- (sum(abs(error1)))/length(y)
ema_1

## [1] 3.666667

Construimos funciones personalizadas para cada indicador.

ema <- function(x){
  l <- length(x)
  s <- sum(abs(x))
  s/l
}

ema(error1)

## [1] 3.666667

ema(error2)

## [1] 3.333333

1.4.2 Error cuadrático medio

Es una medida de precisión que indica que tan dis perso está el modelo de los datos reales pero cuenta con la característica de brindarle mayor peso a los errores más grandes

\[ ECM= \frac{1}{T}\Sigma_{t=1}^T(e_t)^2 \]

ecm <- function(x){
  l <- length(x)
  s <- sum(x^2)
  s/l
}

ecm(error1)

## [1] 21.88889

ecm(error2)

## [1] 20.88889

1.4.3 Error medio absoluto porcentual

Es la medida relativa del EMA

\[ EMAP= \frac{1}{T}\Sigma_{t=1}^T\frac{|e_t|}{y_t}\times100 \]

emap <- function(x,y){
  l <- length(x)
  s <- sum(abs(x)/y)*100
  s/l
}
emap(error1,y)

## [1] 9.856951

emap(error2,y)

## [1] 8.889399

1.4.4 Raíz del error cuadrático medio porcentual

\[ ECMP= \sqrt{\frac{1}{T}\sum_{t=1}^T \left(\frac{e_t}{y_t} \right)^2}\times 100 \]

ecmp <- function(x,y){
  l <- length(x)
  s <- (sum(x/y)^2)
  ss <- sqrt(s/l)*100
  ss
}
ecmp(error1,y)

## [1] 22.36365

ecmp(error2,y)

## [1] 9.784647

1.4.5 U de Theil

\[U = \sqrt{\frac{\frac{1}{T}\sum_{t=1}^Te_t^2}{\frac{1}{T}\sum_{t=1}^Ty_t^2}} \]

UTheil <- function(x,y){
  s <- (sum(x)^2)/(sum(y)^2)
  ss <- sqrt(s)
  ss
}
UTheil(error1,y)

## [1] 0.08012821

UTheil(error2,y)

## [1] 0.03846154

1.4.6 Función agrupada

ts_precision <- function(x, y, 
                         type = c("ecm", "ecmp", "ema", "emap", 
                                  "UTheil"),
                         na.rm = TRUE)
{
  switch(match.arg(type),
  ecm = ecm(x),
  ema = ema(x),
  emap = emap(x, y),
  ecmp = ecmp(x, y),
  UTheil=UTheil(x,y))
}
ema <- function(x){
  l <- length(x)
  s <- sum(abs(x))
  s/l
}
ecm <- function(x){
  l <- length(x)
  s <- sum(x^2)
  s/l
}
emap <- function(x,y){
  l <- length(x)
  s <- sum(abs(x)/y)*100
  s/l
}
ecmp <- function(x,y){
  l <- length(x)
  s <- (sum(x/y)^2)
  ss <- sqrt(s/l)*100
  ss
}
UTheil <- function(x,y){
  s <- (sum(x)^2)/(sum(y)^2)
  ss <- sqrt(s)
  ss
}
ts_precision(error1,y,type="ecm")

## [1] 21.88889

1.4.7 Función conjunta

Construimos una función que además de calcular todos los indicadores de precisión, los pone en una tabla comparativa.

ts_precision_all <- function(x,y){
  ema <- ts_precision(x,y,type="ema")
  emap<- ts_precision(x,y,type="emap")
  ecm <- ts_precision(x,y,type="ecm")
  ecmp <- ts_precision(x,y,type="ecmp")
  u <- ts_precision(x,y,type="UTheil")
  tabla <- as.matrix(data.frame(ema, emap, ecm, ecmp, u))
  colnames(tabla)<- c("EMA", "EMAP", "ECM", "ECMP", "UTheil")
  knitr::kable(
  tabla, booktabs = TRUE,
  caption = "Indicadores de Precisión"
)
  
}

ts_precision_all(error1, y)

Indicadores de Precisión

EMA	EMAP	ECM	ECMP	UTheil
3.666667	9.856951	21.88889	22.36365	0.0801282

ts_precision_all(error2, y)

Indicadores de Precisión

EMA	EMAP	ECM	ECMP	UTheil
3.333333	8.889399	20.88889	9.784647	0.0384615

2 Modelos de Patrones de Comportamiento Subyacentes

Los modelos de series de tiempo de Patrones de Comportamiento Subyacente pretenden recoger las regularidades de la serie capturadas en el comportamiento de una variable a lo largo del tiempo. Descansab en la idea de que los componentes subyacentes se pueden reproducir y nos permiten analizar los datos observados. La metodología es muy básica: se observa un determinado fenóomeno, se analizan sus regularidades o patrones y se construyen modelos basados en dichas regularidades observadas. Descansan en la idea de que toda serie de tiempo se puede descomponer en diferentes elementos subyacentes que forman parte de la serie, de manera indirecta, y que pueden explicar su evolución, a lo largo del tiempo. Los Modelos del Componentes Subyacentes (MCS) no pretenden representar el proceso generador de datos sino, mas bien, describir, de manera muy general, las características fundamentales de las series.

La construcción de que una serie de tiempo se puede construir por elementos superpuestos no es nueva, ya en la edad media se preguntaban por que había fenómenos qe se repetían de manera regular en el tiempo. El punto de partida para la construcción de este tipo de modelos es, por lo tanto, saber identificar la existencia de dichos patrones o regularidades. Este enfoque pretende clasificar los tipos de movimientos que caracterizan una serie de tiempo como tendencia, estacionalidad, ciclo e irregular. Estos cuatro componentes son considerados como los patrones de comportamiento subyacente.

El objetivo del análisis de series de tiempo es dentificar los patrones subyacentes expresados en un modelo uniecuacional.

\[Y_t=f(T_t, C_t, E_t, I_t) \] El componente irregular \(I_t\) representa todos los movimientos no sistmáticos de la serie a lo largo del tiempo. Son las perturbaciones aleatorias que se presentan en la vida cotidiana. Son impredecibles y en suma su efecto se anula entre ellas. La estacionalidad es un patrón que se reproduce al interior de un año y refleja los impulsos del fenómeno asociados con efectos observados en diversas estaciones del año. El ciclo es un comportamiento de mediano plazo que, normalmente son obserbables en series con temporalidades superiores a un año.Por último, la tendencia refleja el comportamiento de largo plazo y representa el movimiento general de la serie una vez que se han eliminado los otros tres efectos.

2.1 Tendencia

La tendencia se define como el comportamiento subyacente de largo plazo en un fenómeno temporal. Normalmente, se realizan ajustes de funciones matemáticas para representar los diferentes tipos de tendencia.

2.1.1 Tendencia lineal

Supongamos que estamos interesados en evaluar el comportamiento del Índice Nacional de Precios al Consumidor (INPC). El INPC se encuentra en el Banco de Información Económica del INEGI¹ y su identificador es 628194.

inegi_id <- "628194"
INPC <- inegi_series(inegi_id, token, database = "BIE")
INPC <- INPC[1:367,]
INPC <- ts_invert(INPC)
INPC <- ts(INPC, frequency = 12, start = c(1994,1))
tail(INPC)

## [1] 133.681 134.065 134.336 134.087 134.594 136.003

Como primer paso, observamos el fenómeno:

ts.plot(INPC)

Como se puede apreciar, hay una línea con pocas variaciones. De manera tal que podemos plantear un modelo matemático que incluya la tendencia lineal \(T_t\) y el componente irregular \(I_t\).

\[ Y_t = f(t_t, I_t)\] Si el modelo fuera aditivo, tomaría la siguiente forma

\[Y_t= T_t+I_t\] En particular, el componente de tendencia puede ser modelado como una ecuación lineal con \(\alpha\) como intercepto y \(\beta\) como pendiente

\[T_t= \alpha + \beta t \] En el contexto de las series de tiempo, el componente irregular se identifica con el término de perturbación estócastica, se llama innovación y su nomenclatura es \(a_t\)

\[Y_t =\alpha + \beta t+a_t \]

Para modelar la tendencia, es necesario construir la variable tiempo

tiempo <- time(INPC)
plot(tiempo)

\[y_t=\alpha+\beta t+\epsilon_t \] Para modelar la tendencia es necesario construir un modelo econométrico para la función lineal propuesta

\[\hat{y_t}=\hat{\alpha}+\hat{\beta} t \]

lm(INPC~tiempo)

## 
## Call:
## lm(formula = INPC ~ tiempo)
## 
## Coefficients:
## (Intercept)       tiempo  
##   -6856.630        3.448

El resumen estadístico del modelo se imprime con la función summary()

mod0 <- lm(INPC~tiempo)
summary(mod0)

## 
## Call:
## lm(formula = INPC ~ tiempo)
## 
## Residuals:
##    Min     1Q Median     3Q    Max 
## -6.868 -3.113 -1.191  2.739 12.314 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -6.857e+03  4.961e+01  -138.2   <2e-16 ***
## tiempo       3.448e+00  2.469e-02   139.6   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 4.176 on 365 degrees of freedom
## Multiple R-squared:  0.9816, Adjusted R-squared:  0.9816 
## F-statistic: 1.95e+04 on 1 and 365 DF,  p-value: < 2.2e-16

A partir del resumen estadístico del modelo de tendencia lineal se pueden almacenar los coeficientes de la regresión.

smod0 <- summary(mod0)
a <- coef(mod0)[[1]]
a

## [1] -6856.63

b <- coef(mod0)[[2]]
b

## [1] 3.447922

Con base en dichos coeficientes, se realiza un ajuste con los datos observados y los datos estimados por medio de los coeficientes del modelo lineal.

plot(INPC)
abline(a,b, col="salmon")

El análisis de series de tiempo tiene dos grandes utilidades, por un lado es factible evaluar el comportamiento histórico del fenómeno en cuestión para analizar el comportamiento a lo largo del tiempo. Por otro lado, con base en la reproducción de los patrones subyacentes, es factible construir, si existe estabilidad paramétrica en los patrones subyacentes, un comportamiento que pemrita reproducir dichos patrones hacia el futuro. Por medio de la librería forecast se puede alcanzar este objetivo.

#install.packages("forecast")
library(forecast)

Construimos la tendencia lineal con los coeficientes y la trasnformamos a serie de tiempo con la temporalidad y la frecuencia de la variable observada y construimos el prónóstico para un horizonte temporal de cinco momentos en el futuro.

Al imprimir el pronóstico se precia el valor estimado del valor futuro, así como sus intervalos con 80% y 90% de confianza.

tendenciaL <- ts(a+b*tiempo, frequency = 12, start = c(1994,1))
pronostico <- forecast(tendenciaL, h=5)
pronostico

##          Point Forecast    Lo 80    Hi 80    Lo 95    Hi 95
## Aug 2024       123.9759 123.9759 123.9759 123.9759 123.9759
## Sep 2024       124.2633 124.2633 124.2633 124.2633 124.2633
## Oct 2024       124.5506 124.5506 124.5506 124.5506 124.5506
## Nov 2024       124.8379 124.8379 124.8379 124.8379 124.8379
## Dec 2024       125.1253 125.1253 125.1253 125.1253 125.1253

Adicionalmente se puede graficar el pronóstico de la tendencia junto con la tendencia histórica

plot(forecast(tendenciaL))

Una manera mucho mas atractiva de mostrar los resultados es por medio de la función autoplot() de la librería forecast y las opciones gráficas de la librería ggplot2.

library(forecast)
library(ggplot2)
inpc <- window(INPC,start=c(1994,1))
fit1 <- tendenciaL
autoplot(inpc) +
autolayer(pronostico, series="Pronóstico")+
autolayer(fit1, series="Tendencia lineal") +
  xlab("Año") +
  ylab("INPC (índice)") +
  ggtitle("Indice Nacional de Precios al Consumidor") +
  guides(colour=guide_legend(title="Pronóstico"))

Del modelo se desprende que, caeteris paribus, el incremento mensual del Índice Nacional de Precios al Consumidor es de 3.4479223 y cuenta con un intercepto de -6856.6300375

2.1.2 Tendencia cuadrática

Modelar una tendencia que , en principio, parece ser cuadrática. La variable es la Tasa de condiciones críticas de ocupación (TCCO) (Porcentaje), cuyo identificador es 444608. La tasa de condiciones críticas de ocupación representa el segmento de la población económicamente activa que cuenta con alguna de las siguientes características:

1. Población que trabaja menos de 35 horas a la semana
2. Población que trabaja mas de 35 horas a la semana y recibe hasta un salario mínimo.
3. Población que trabaja mas de 48 horas y hasta menos de dos salarios mínimos.

inegi_id <- "444608"
TCCO <- inegi_series(inegi_id, token, database = "BIE")
TCCO <- ts_invert(TCCO)
TCCO <- ts(TCCO, frequency = 12, start = c(2005,1))
ts.plot(TCCO)

\[y_t=\beta_0+\beta_1 t+ \beta_3 t^2 +\epsilon_t\]

tiempo <- time(TCCO)
tiempo2 <- tiempo^2
mod1 <- tslm(TCCO~tiempo+tiempo2)
mod1a<- tslm(TCCO~trend+I(trend^2))
b0<-coef(mod1)[[1]]
b1<-coef(mod1)[[2]]
b2<-coef(mod1)[[3]]
tendenciaC <- b0+b1*tiempo+b2*tiempo2
plot(tendenciaC, col="aquamarine2")

ts.plot(TCCO, tendenciaC, col=1:2)

plot(forecast(mod1a, h=12))

2.1.3 Tendencia Cúbica

Tasa de ocupación parcial y desocupación 1 (TOPD1) (Porcentaje), con identificador 444604. La tasa de ocupación parcial y desocupación es el porcentaje de la población económicamente activa que se encuentra desocupada, más la ocupada que trabajó menos de 15 horas en la semana.

inegi_id <- "444604"
TOPD1 <- inegi_series(inegi_id, token, database = "BIE")
TOPD1 <- ts_invert(TOPD1)
TOPD1 <- ts(TOPD1, frequency = 12, start = c(2005,1))
ts.plot(TOPD1)

tiempo <- seq(1, length(TOPD1), 1)
tiempo2 <- tiempo^2
tiempo3 <- tiempo^3
mod3 <- lm(TOPD1~tiempo+tiempo2+tiempo3)
mod3a <- tslm(TOPD1~trend++I(trend^2)++I(trend^3))
summary(mod3)

## 
## Call:
## lm(formula = TOPD1 ~ tiempo + tiempo2 + tiempo3)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -1.7320 -0.6310 -0.1026  0.3904  5.4030 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  8.481e+00  2.519e-01  33.662  < 2e-16 ***
## tiempo       7.367e-02  9.110e-03   8.087 3.30e-14 ***
## tiempo2     -6.169e-04  8.847e-05  -6.973 3.14e-11 ***
## tiempo3      1.363e-06  2.434e-07   5.600 5.98e-08 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.9564 on 234 degrees of freedom
## Multiple R-squared:  0.3492, Adjusted R-squared:  0.3409 
## F-statistic: 41.85 on 3 and 234 DF,  p-value: < 2.2e-16

b0 <- coef(mod3)[[1]]
b1 <- coef(mod3)[[2]]
b2 <- coef(mod3)[[3]]
b3 <- coef(mod3)[[4]]
tendencia <- b0+b1*tiempo+b2*tiempo2+b3*tiempo3
ts.plot(TOPD1, tendencia, col=1:2)

forecast(tendencia, h=5)

##     Point Forecast    Lo 80    Hi 80    Lo 95    Hi 95
## 239       9.453163 9.452620 9.453705 9.452333 9.453992
## 240       9.463922 9.462744 9.465100 9.462120 9.465724
## 241       9.474466 9.472527 9.476405 9.471500 9.477432
## 242       9.484800 9.481996 9.487603 9.480512 9.489087
## 243       9.494926 9.491170 9.498682 9.489182 9.500671

plot(forecast(tendencia, h=12))

plot(forecast(mod3a, h=12))

2.1.4 Tendencia Polinomial

Tasa de presión general (TPRG) (Porcentaje), 444605 polinomial

inegi_id <- "444605"
TPRG <- inegi_series(inegi_id, token, database = "BIE")
TPRG <- ts_invert(TPRG)
TPRG <- ts(TPRG, frequency = 12, start = c(2005,1))
plot(TPRG, type="l")

tiempo <- seq(1, length(TPRG), 1)
tiempo2 <- tiempo^2
tiempo3 <- tiempo^3
tiempo4 <- tiempo^4
mod4 <- lm(TPRG~tiempo+tiempo2+tiempo3+tiempo4)
mod4a <- tslm(TPRG~trend+I(trend^2)+I(trend^3)+I(trend^4))
summary(mod4)

## 
## Call:
## lm(formula = TPRG ~ tiempo + tiempo2 + tiempo3 + tiempo4)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -1.69922 -0.72949 -0.03688  0.59556  2.43773 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  6.725e+00  2.864e-01  23.485  < 2e-16 ***
## tiempo       5.863e-02  1.653e-02   3.548  0.00047 ***
## tiempo2     -5.334e-04  2.804e-04  -1.903  0.05833 .  
## tiempo3      1.343e-06  1.761e-06   0.762  0.44659    
## tiempo4     -7.905e-10  3.655e-09  -0.216  0.82896    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.8613 on 233 degrees of freedom
## Multiple R-squared:  0.4921, Adjusted R-squared:  0.4834 
## F-statistic: 56.44 on 4 and 233 DF,  p-value: < 2.2e-16

b0 <- coef(mod4)[[1]]
b1 <- coef(mod4)[[2]]
b2 <- coef(mod4)[[3]]
b3 <- coef(mod4)[[4]]
b4 <- coef(mod4)[[5]]
tendencia <- b0+b1*tiempo+b2*tiempo2+b3*tiempo3+b4*tiempo4
ts.plot(TPRG, tendencia, col=1:2)

plot(forecast(mod4a, h=12))

2.1.5 Tendencia exponencial

Finanzas públicas > Ingresos y egresos brutos por entidad federativa > Nacional 31 entidades federativas (no incluye Ciudad de México) > Ingresos > Impuestos Clave:441566

inegi_id <- "441566"
impuestos <- inegi_series(inegi_id, token, database = "BIE")
impuestos <- ts_invert(impuestos)
impuestos <- ts(impuestos, start = 1989)
ts.plot(impuestos )

\[y_t=T_tI_t \] \[y_t=e^{\beta_0+\beta_1t}e^{a_t} \] \[y_t=e^{\beta_0+\beta_1t+a_t} \] \[ln(y_t)= \beta_0+\beta_1t+a_t \]

limpuestos<- log(impuestos)
tiempo <- time(impuestos)
mod5 <- lm(limpuestos~tiempo)
mod5a <- tslm(log(impuestos)~trend)

b0 <- coef(mod5)[[1]]
b1 <- coef(mod5)[[2]]
sub_tendencia <- b0+b1*tiempo
tendencia <- exp(sub_tendencia)
ts.plot(impuestos, tendencia, col=1:2)

plot(forecast(mod5a, h=4))

2.1.6 Curvas de crecimiento

## Rows: 364 Columns: 5
## ── Column specification ────────────────────────────────────────────────────────
## Delimiter: ","
## chr (1): fecha
## dbl (4): diarios, diarios_cdmx, totales, totales_cdmx
## 
## ℹ Use `spec()` to retrieve the full column specification for this data.
## ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.

dias <- seq(1,length(mexicoc$totales),1)
totales <- mexicoc$totales
datos <-data.frame(dias, totales)
dias2 <- seq(1,387,1)
totales2 <- mexicoc$totales[1:387]
datos2 <-data.frame(dias2, totales2)

ajuste <-nls(totales~a*exp(-exp(-((dias-b)/c))), 
             start=list(a=500000, b=100, c=4), trace = T)#      

## 7.646194e+13 (6.77e-01): par = (5e+05 100 4)
## 4.346055e+13 (1.01e+00): par = (798186.9 106.3715 22.29686)
## 2.116812e+13 (2.98e+00): par = (1029113 135.3097 77.76)
## 1.444033e+13 (4.24e+00): par = (1812924 245.4666 158.3111)
## 9.841359e+12 (3.32e+00): par = (3149755 333.1296 170.8134)
## 8.787634e+11 (2.82e-01): par = (3653680 320.6704 175.9763)
## 8.210666e+11 (4.68e-02): par = (3248955 302.4053 159.7596)
## 8.193896e+11 (7.27e-03): par = (3339346 306.6954 163.1357)
## 8.193508e+11 (1.16e-03): par = (3336786 306.5205 162.8819)
## 8.193499e+11 (2.11e-04): par = (3339425 306.6546 162.992)
## 8.193499e+11 (3.89e-05): par = (3339163 306.641 162.9785)
## 8.193499e+11 (7.22e-06): par = (3339238 306.6448 162.9818)

sajuste<- summary(ajuste)
a<- sajuste$coefficients[1]
b<- sajuste$coefficients[2]
c<- sajuste$coefficients[3]
c(a,b,c) 

## [1] 3339237.8518     306.6448     162.9818

f<-seq(1,1000,1)
contagiados<-a*exp(-exp(-((f-b)/c)))
#format(contagiados, scientific=F) 
options(scipen=999)
contagiados <- ceiling(contagiados)
plot(f,contagiados, type = "h", 
     xlab = "Días", ylab = "Contagios acumulados", col="blue")

dias2 <- f
contagios <- data.frame(dias2, contagiados)
tail(contagios)

##      dias2 contagiados
## 995    995     3290684
## 996    996     3290979
## 997    997     3291272
## 998    998     3291564
## 999    999     3291853
## 1000  1000     3292141

g <- data.frame(f, contagiados)
dias3 <- d<-seq(1,length(mexicoc),1)
totales3 <- c(contagiados[1:length(mexicoc)])
datos3 <-data.frame(dias3, totales3)

library(ggplot2)
ggplot(g, aes(x = f, y = contagiados, colour = "Curva Gompertz") ) +
   geom_line()+ theme(legend.position = "bottom")+
   geom_point(data = datos, aes(x=dias, y = totales), size = 1, color="orange")+ theme(legend.position = "bottom")+
   #geom_line(data=deriv, aes(x=f1, y=derivada, color="brown"))+ 
   xlab("Días desde el pacientel 0") + ylab("Número de contagios confirmados")+
   ggtitle("Modelo de predicción de contagios confirmados de Covid-19.")+
   labs(subtitle = "Ajuste de una función de crecimiento tipo Gompertz",
   caption = "Elaborado por el Laboratorio de Desigualdad Socioeconómica Regional DeSERLab, UAM-I.")+ 
   geom_hline(yintercept=3300000, color="red", show.legend =TRUE, linetype="dashed")+
   geom_vline(xintercept=1000, color="red", show.legend=TRUE, linetype="dashed")+
   annotate("text", x=250, y=3100000, label="Valor asintótico: 10,965,030")+
   annotate("text", x=600, y=2000000, label="2024-02-19")+
   annotate("text", x=650, y=2200000, label="Fecha probable del Plateau")

2.2 Estacionalidad

2.2.1 Estacionalidad aditiva

Comportamiento que se presenta de manera regular y sistemática

737124 112 Cría y explotación de animales (Índice base 2018=100) Indicadores económicos de coyuntura > Indicador global de la actividad económica, base 2018 > Series originales > Índice de volumen físico > Actividades primarias estacionalidad aditiva

inegi_id <- "737124"
IVFEA <- inegi_series(inegi_id, token, database = "BIE")
IVFEA <- ts_invert(IVFEA)
IVFEA <- ts(IVFEA, frequency = 12, start = c(1993,1))
ts.plot(IVFEA)

Seasonal Plot:

library(ggplot2)
ggseasonplot(IVFEA, year.labels=TRUE, year.labels.left=TRUE) +
  ylab("Índice base 2018=100") +
  ggtitle("Seasonal plot: Cría y explotación de animales")+
  xlab("Mes")

Polar Seasonal Plot:

ggseasonplot(IVFEA, polar=TRUE) +
  ylab("Índice base 2018=100") +
  ggtitle("Polar seasonal plot: Cría y explotación de animales")+
  xlab("Mes")

Subseries Plot:

ggsubseriesplot(IVFEA) +
  ylab("Índice base 2018=100") +
  ggtitle("Subseries seasonal plot: Cría y explotación de animales")+
  xlab("Mes")

Lag plots

ivfea <- window(IVFEA, start=1993)
gglagplot(ivfea)

Here the colours indicate the quarter of the variable on the vertical axis. The lines connect points in chronological order. The relationship is strongly positive at lags 4 and 8, reflecting the strong seasonality in the data. The negative relationship seen for lags 2 and 6 occurs because peaks (in Q4) are plotted against troughs (in Q2)

The window() function used here is very useful when extracting a portion of a time series. In this case, we have extracted the data from ausbeer, beginning in 1992.

library(forecast)
dummies <- seasonaldummy(IVFEA)
head(dummies)

##      Jan Feb Mar Apr May Jun Jul Aug Sep Oct Nov
## [1,]   1   0   0   0   0   0   0   0   0   0   0
## [2,]   0   1   0   0   0   0   0   0   0   0   0
## [3,]   0   0   1   0   0   0   0   0   0   0   0
## [4,]   0   0   0   1   0   0   0   0   0   0   0
## [5,]   0   0   0   0   1   0   0   0   0   0   0
## [6,]   0   0   0   0   0   1   0   0   0   0   0

m6 <- lm(IVFEA~dummies)
plot(m6$fitted.values, type="l")

ts.plot(IVFEA, m6$fitted.values, col=1:2)

m6a <- tslm(IVFEA~ trend + season)
ts.plot(IVFEA, m6a$fitted.values, col=1:2)

plot(forecast(m6a, h=12))

2.2.2 Estacionalidad multiplicativa

737123 111 Agricultura (Índice base 2018=100) Indicadores económicos de coyuntura > Indicador global de la actividad económica, base 2018 > Series originales > Índice de volumen físico > Actividades primarias estacionalidad multiplicativa

inegi_id <- "737123"
IVFAgr <- inegi_series(inegi_id, token, database = "BIE")
IVFAgr <- ts_invert(IVFAgr)
IVFAgr <- ts(IVFAgr, frequency = 12, start = c(1993,1))
ts.plot(IVFAgr)

tiempo<- time(IVFAgr)
dummies <- seasonaldummy(IVFAgr)
mod7 <- lm(IVFAgr~tiempo+dummies)
summary(mod7)

## 
## Call:
## lm(formula = IVFAgr ~ tiempo + dummies)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -24.948  -5.867  -0.866   5.515  33.595 
## 
## Coefficients:
##                Estimate  Std. Error t value             Pr(>|t|)    
## (Intercept) -3186.91639   105.41709 -30.231 < 0.0000000000000002 ***
## tiempo          1.63915     0.05247  31.241 < 0.0000000000000002 ***
## dummiesJan    -19.03954     2.36509  -8.050   0.0000000000000116 ***
## dummiesFeb    -34.21718     2.36505 -14.468 < 0.0000000000000002 ***
## dummiesMar    -42.03273     2.36502 -17.773 < 0.0000000000000002 ***
## dummiesApr    -33.09039     2.36500 -13.992 < 0.0000000000000002 ***
## dummiesMay    -18.02952     2.36499  -7.624   0.0000000000002122 ***
## dummiesJun    -13.54341     2.36499  -5.727   0.0000000213018423 ***
## dummiesJul    -24.23038     2.36499 -10.245 < 0.0000000000000002 ***
## dummiesAug    -47.34569     2.36500 -20.019 < 0.0000000000000002 ***
## dummiesSep    -54.78588     2.36502 -23.165 < 0.0000000000000002 ***
## dummiesOct    -43.00850     2.38370 -18.043 < 0.0000000000000002 ***
## dummiesNov     -2.80241     2.38369  -1.176                 0.24    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 9.385 on 368 degrees of freedom
## Multiple R-squared:  0.8558, Adjusted R-squared:  0.8511 
## F-statistic:   182 on 12 and 368 DF,  p-value: < 0.00000000000000022

ts.plot(IVFAgr, mod7$fitted.values, col=1:2)

mod7a <- lm(log(IVFAgr)~tiempo+dummies)
summary(mod7a)

## 
## Call:
## lm(formula = log(IVFAgr) ~ tiempo + dummies)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.31782 -0.05770 -0.00164  0.06746  0.35158 
## 
## Coefficients:
##                Estimate  Std. Error t value             Pr(>|t|)    
## (Intercept) -36.6805133   1.1257796 -32.582 < 0.0000000000000002 ***
## tiempo        0.0205659   0.0005603  36.704 < 0.0000000000000002 ***
## dummiesJan   -0.1968438   0.0252575  -7.793   0.0000000000000674 ***
## dummiesFeb   -0.3788068   0.0252571 -14.998 < 0.0000000000000002 ***
## dummiesMar   -0.4855111   0.0252568 -19.223 < 0.0000000000000002 ***
## dummiesApr   -0.3677764   0.0252565 -14.562 < 0.0000000000000002 ***
## dummiesMay   -0.1877426   0.0252564  -7.433   0.0000000000007487 ***
## dummiesJun   -0.1308223   0.0252564  -5.180   0.0000003666968376 ***
## dummiesJul   -0.2508066   0.0252564  -9.930 < 0.0000000000000002 ***
## dummiesAug   -0.5766381   0.0252565 -22.831 < 0.0000000000000002 ***
## dummiesSep   -0.7130341   0.0252568 -28.231 < 0.0000000000000002 ***
## dummiesOct   -0.5170141   0.0254562 -20.310 < 0.0000000000000002 ***
## dummiesNov   -0.0352248   0.0254561  -1.384                0.167    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.1002 on 368 degrees of freedom
## Multiple R-squared:  0.8926, Adjusted R-squared:  0.8891 
## F-statistic: 254.9 on 12 and 368 DF,  p-value: < 0.00000000000000022

ts.plot(IVFAgr, exp(mod7a$fitted.values), col=1:2)

lambda= Box-Cox transformation parameter. If lambda=“auto”, then a transformation is automatically selected using BoxCox.lambda. The transformation is ignored if NULL. Otherwise, data transformed before model is estimated.

mod7b <- tslm(IVFAgr~trend+season, lambda = -0.6)
ts.plot(IVFAgr,mod7b$fitted.values, col=1:2)

plot(forecast(mod7b, h=12))

2.3 Ciclo

Moviemiento de mediano plazo

737126 21 Minería (Índice base 2018=100) Indicadores económicos de coyuntura > Indicador global de la actividad económica, base 2018 > Series originales > Índice de volumen físico > Actividades secundarias

inegi_id <- "737126"
IVFMin <- inegi_series(inegi_id, token, database = "BIE")
IVFMin <- ts_invert(IVFMin)
IVFMin <- ts(IVFMin, frequency = 12, start = c(1993,1))
ts.plot(IVFMin)

head(cycle(IVFMin),n=12)

##      Jan Feb Mar Apr May Jun Jul Aug Sep Oct Nov Dec
## 1993   1   2   3   4   5   6   7   8   9  10  11  12

aggregate.ts(IVFMin, FUN=mean)

## Time Series:
## Start = 1993 
## End = 2023 
## Frequency = 1 
##  [1] 106.77938 106.63936 104.92355 116.46393 123.60532 126.04125 124.00032
##  [8] 127.91148 130.76937 129.10725 133.86573 136.19819 136.93546 135.83764
## [15] 134.25355 128.09959 123.31333 125.66709 126.42954 128.72318 128.72786
## [22] 126.62119 120.72698 116.26661 106.90804 100.00000  94.62836  94.22457
## [29]  96.64103 100.60956 100.74670

plot(aggregate(IVFMin, FUN=mean))

boxplot(IVFMin~cycle(IVFMin))

plot(forecast(aggregate(IVFMin, FUN=mean),h=4))

3 Descomposiciones, Filtros y Suavizamientos

3.1 Descomponer la serie en sus cuatro componentes

778083 -778114 Aguascalientes (Índice de volumen físico 2018 = 100) Indicadores económicos de coyuntura > Indicador trimestral de la actividad económica estatal (ITAEE), base 2018 > Series Originales > Actividades primarias > 111-112 - Agricultura. Cría y explotación de animales > Índice

inegi_id <- "778083"
IVFEAgs <- inegi_series(inegi_id, token, database = "BIE")
IVFEAgs <- ts_invert(IVFEAgs)
IVFEAgs <- ts(IVFEAgs, frequency = 4, start = c(2003,1))
ts.plot(IVFEAgs)

descomposicion <- decompose(IVFEAgs)
head(descomposicion$seasonal)

##            Qtr1       Qtr2       Qtr3       Qtr4
## 2003 -13.799073  -3.184060  11.267782   5.715351
## 2004 -13.799073  -3.184060

head(descomposicion$trend)

##          Qtr1     Qtr2     Qtr3     Qtr4
## 2003       NA       NA 55.69887 55.74841
## 2004 55.22556 54.89617

head(descomposicion$random)

##           Qtr1      Qtr2      Qtr3      Qtr4
## 2003        NA        NA  0.786847  1.302763
## 2004  2.177347 -3.870247

head(descomposicion$figure)

## [1] -13.799073  -3.184060  11.267782   5.715351

plot(descomposicion)

dummy <- descomposicion$seasonal
tendencia <- descomposicion$trend
y <- tendencia +dummy
ts.plot(IVFEAgs,y, col=1:2)

descomposicion2 <- data.frame(cbind(l=log(IVFEAgs),
                                    t=time(log(IVFEAgs)),
                                    c=cycle(log(IVFEAgs)),
                                    s=factor(descomposicion$seasonal)))
mod8 <- lm(l~t+c+s, data=descomposicion2)
smod8 <- summary(mod8)
yhat <- mod8$fitted.values
ts.plot(IVFEAgs, exp(yhat), col=1:2)

3.2 Promedios Móviles

SMA= Simple Moving Average EMA= Exponential Moving Average WMA = Weigthed Moving Average DEMA = Differential Exponential Moving Average EVWMA = Elastic Volume Weigthed Moving Average ZLEMA = Zero Lag Exponential Moving Avergae WVMA = Volume Weigthed Moving Average HMA = Diferencia de dos WMA ALMA = Filtro GAussiano

Precio de Apertura Precio de Cierre Precio Máximo Precio Mínimo Volumen

library(TSA)

## Registered S3 methods overwritten by 'TSA':
##   method       from    
##   fitted.Arima forecast
##   plot.Arima   forecast

## 
## Attaching package: 'TSA'

## The following object is masked from 'package:readr':
## 
##     spec

## The following objects are masked from 'package:stats':
## 
##     acf, arima

## The following object is masked from 'package:utils':
## 
##     tar

#install.packages("TTR")
library(TTR)
data(airpass)
ts.plot(IVFEAgs, SMA(IVFEAgs, n=10), col=1:2)

tail(SMA(IVFEAgs, n=10))

##          Qtr1     Qtr2     Qtr3     Qtr4
## 2023 104.6562 105.1373 109.6169 109.2868
## 2024 105.3327 105.2148

ts.plot(IVFEAgs, SMA(IVFEAgs, n=10),
        SMA(IVFEAgs, n=3),
        SMA(IVFEAgs, n=25),
        SMA(IVFEAgs, n=40),col=1:5)

library(forecast)
plot(forecast(SMA(IVFEAgs, n=5)))

## Warning in ets(object, lambda = lambda, biasadj = biasadj,
## allow.multiplicative.trend = allow.multiplicative.trend, : Missing values
## encountered. Using longest contiguous portion of time series

3.3 Filtros

Official statistics agencies (such as the US Census Bureau and the Australian Bureau of Statistics) are responsible for a large number of official economic and social time series. These agencies have developed their own decomposition procedures which are used for seasonal adjustment. Most of them use variants of the X-11 method, or the SEATS method, or a combination of the two. These methods are designed specifically to work with quarterly and monthly data, which are the most common series handled by official statistics agencies. They will not handle seasonality of other kinds, such as daily data, or hourly data, or weekly data. We will use the latest implementation of this group of methods known as “X-13ARIMA-SEATS”. For the methods discussed in this section, you will need to have installed the seasonal package in R.

X-11 method The X-11 method originated in the US Census Bureau and was further developed by Statistics Canada. It is based on classical decomposition, but includes many extra steps and features in order to overcome the drawbacks of classical decomposition that were discussed in the previous section. In particular, trend-cycle estimates are available for all observations including the end points, and the seasonal component is allowed to vary slowly over time. X-11 also handles trading day variation, holiday effects and the effects of known predictors. There are methods for both additive and multiplicative decomposition. The process is entirely automatic and tends to be highly robust to outliers and level shifts in the time series. The details of the X-11 method are described in Dagum & Bianconcini (2016).

plot(IVFEAgs)

library(forecast)
lambda <- BoxCox.lambda(IVFEAgs)
IVFEAgs.fit <- ar(BoxCox(IVFEAgs,lambda))
plot(forecast(IVFEAgs,h=20,lambda=lambda))

fit <- ets(IVFEAgs, model = "ZZZ", allow.multiplicative.trend = TRUE)
plot(forecast(fit))

library(ggplot2)

# Using Fourier series for a "ts" object
# K is chosen to minimize the AICc
IVFEAgs.model  <- auto.arima(IVFEAgs, xreg=fourier(IVFEAgs,K=2), seasonal=FALSE)
IVFEAgs.fcast <- forecast(IVFEAgs.model, xreg=fourier(IVFEAgs, K=2, h=8))
autoplot(IVFEAgs.fcast) + xlab("Year")

IVFEAgs.fcast <- rwf(IVFEAgs, h=4)
plot(IVFEAgs.fcast)

plot(naive(IVFEAgs,h=8),include=200)

plot(snaive(IVFEAgs))

plot(ses(IVFEAgs))

fcast <- holt(IVFEAgs)
plot(fcast)

IVFEAgs.fcast <- hw(IVFEAgs,h=8)
plot(IVFEAgs.fcast)

fcast <- ets(IVFEAgs)
plot(fcast)

HoltWinters(IVFEAgs)

## Holt-Winters exponential smoothing with trend and additive seasonal component.
## 
## Call:
## HoltWinters(x = IVFEAgs)
## 
## Smoothing parameters:
##  alpha: 0.05804195
##  beta : 0.5067056
##  gamma: 0.446458
## 
## Coefficients:
##            [,1]
## a  106.90612196
## b    0.07312874
## s1  18.50072372
## s2   1.44704963
## s3 -22.30953812
## s4   7.59852503

plot(HoltWinters(IVFEAgs))

IVFEAgsHW <- HoltWinters(IVFEAgs)
forecast(IVFEAgsHW)

##         Point Forecast     Lo 80     Hi 80     Lo 95     Hi 95
## 2024 Q3      125.47997 118.18408 132.77587 114.32186 136.63809
## 2024 Q4      108.49943 101.17569 115.82317  97.29873 119.70013
## 2025 Q1       84.81597  77.44276  92.18918  73.53962  96.09232
## 2025 Q2      114.79716 107.34712 122.24720 103.40331 126.19101
## 2025 Q3      125.77249 117.14546 134.39952 112.57858 138.96640
## 2025 Q4      108.79194 100.03611 117.54778  95.40104 122.18284
## 2026 Q1       85.10849  76.18705  94.02992  71.46433  98.75264
## 2026 Q2      115.08968 105.96283 124.21653 101.13136 129.04799

plot(forecast(IVFEAgsHW))

library(quantmod)
getFX("USD/MXN")

## [1] "USD/MXN"

plot(HoltWinters(USDMXN[,1], gamma = FALSE))

hh <- ts(USDMXN)
hhh<- ts(hh[,1])
plot(HoltWinters(hhh, gamma = FALSE))

plot(forecast(HoltWinters(hhh, gamma = FALSE)))

3.4 Herramientas adicionales de pronóstico

702255 Total industrias manufactureras (Índice base 2018=100) Indicadores económicos de coyuntura > Encuesta mensual de la industria manufacturera (EMIM). Serie 2018 > Series Originales > Índices de remuneraciones medias reales por persona > Índice > Índice de salarios medios reales por obrero

inegi_id <- "702255"
ISMRO <- inegi_series(inegi_id, token, database = "BIE")
ISMRO <- ts_invert(ISMRO)
ISMRO <- ts(ISMRO, frequency = 12, start = c(2007,1))
ts.plot(ISMRO)

ISMRO <- window(ISMRO,start=2007)
ISMRO1 <- meanf(ISMRO,h=10)
ISMRO2 <- rwf(ISMRO,h=10)
ISMRO3 <- snaive(ISMRO,h=10)
ISMRO4 <- croston(ISMRO,h=10)
ISMRO5 <- stlf(ISMRO,h=10)
ISMRO6 <- ses(ISMRO,h=10)
ISMRO7 <- holt(ISMRO,h=10)
ISMRO8 <- hw(ISMRO,h=10)
ISMRO9 <- splinef(ISMRO,h=10)
ISMRO10 <- thetaf(ISMRO,h=10)
autoplot(window(ISMRO, start=2007)) +
  autolayer(ISMRO1, series="Mean", PI=FALSE) +
  xlab("Mes") + ylab("Índice base 2018=100") +
  ggtitle("pronóstico de Índice de salarios medios reales por obrero") +
  guides(colour=guide_legend(title="Pronóstico"))

autoplot(window(ISMRO, start=2007)) +
  autolayer(ISMRO2, series="Naïve", PI=FALSE) +
  xlab("Mes") + ylab("Índice base 2018=100") +
  ggtitle("pronóstico de Índice de salarios medios reales por obrero") +
  guides(colour=guide_legend(title="Pronóstico"))

autoplot(window(ISMRO, start=2007)) +
  autolayer(ISMRO3, series="Seasonal naïve", PI=FALSE) +
  xlab("Mes") + ylab("Índice base 2018=100") +
  ggtitle("pronóstico de Índice de salarios medios reales por obrero") +
  guides(colour=guide_legend(title="Pronóstico"))

autoplot(window(ISMRO, start=2007)) +
  autolayer(ISMRO4, series="Croston´ Method", PI=FALSE) +
  xlab("Mes") + ylab("Índice base 2018=100") +
  ggtitle("pronóstico de Índice de salarios medios reales por obrero") +
  guides(colour=guide_legend(title="Pronóstico"))

autoplot(window(ISMRO, start=2007)) +
  autolayer(ISMRO5, series="Loess descomposition", PI=FALSE) +
  xlab("Mes") + ylab("Índice base 2018=100") +
  ggtitle("pronóstico de Índice de salarios medios reales por obrero") +
  guides(colour=guide_legend(title="Pronóstico"))

autoplot(window(ISMRO, start=2007)) +
  autolayer(ISMRO6, series="Exponencial Smoothing", PI=FALSE) + 
  xlab("Mes") + ylab("Índice base 2018=100") +
  ggtitle("pronóstico de Índice de salarios medios reales por obrero") +
  guides(colour=guide_legend(title="Pronóstico"))

autoplot(window(ISMRO, start=2007)) +
  autolayer(ISMRO7, series="Holt", PI=FALSE) + 
  xlab("Mes") + ylab("Índice base 2018=100") +
  ggtitle("pronóstico de Índice de salarios medios reales por obrero") +
  guides(colour=guide_legend(title="Pronóstico"))

autoplot(window(ISMRO, start=2007)) +
  autolayer(ISMRO8, series="Holt-Winters", PI=FALSE) + 
  xlab("Mes") + ylab("Índice base 2018=100") +
  ggtitle("pronóstico de Índice de salarios medios reales por obrero") +
  guides(colour=guide_legend(title="Pronóstico"))

autoplot(window(ISMRO, start=2007)) +
  autolayer(ISMRO9, series="cubic SPLine", PI=FALSE) + 
  xlab("Mes") + ylab("Índice base 2018=100") +
  ggtitle("pronóstico de Índice de salarios medios reales por obrero") +
  guides(colour=guide_legend(title="Pronóstico"))

autoplot(window(ISMRO, start=2007)) +
  autolayer(ISMRO10, series="Theta Method", PI=FALSE) +
  xlab("Mes") + ylab("Índice base 2018=100") +
  ggtitle("pronóstico de Índice de salarios medios reales por obrero") +
  guides(colour=guide_legend(title="Pronóstico"))

702140 - 702160 311 Industria alimentaria (Índice base 2018=100) Indicadores económicos de coyuntura > Encuesta mensual de la industria manufacturera (EMIM). Serie 2018 > Series Originales > Índice de personal ocupado > Índice > Índice de personal ocupado total

inegi_id <- "702140"
IPOAlim <- inegi_series(inegi_id, token, geography = "00", database = "BIE")
IPOAlim <- ts_invert(IPOAlim)
IPOAlim <- ts(IPOAlim, frequency = 12, start = c(2007,1))
ts.plot(IPOAlim)

897 Total (Millones de dólares) Indicadores económicos de coyuntura > Balanza comercial de mercancías de México > Series originales > Saldo estacionaria

4 Modelos SARIMA

4.1 Fenómenos Estacionarios

4.1.1 Modelos AR

Indicadores económicos de coyuntura > Tasas de ocupación, desocupación y subocupación (resultados mensuales de la ENOE, 15 años y más) > Series originales > Población total > Población ocupada > Por posición en la ocupación > Trabajadores por cuenta propia Clave:444568

inegi_id <- "444568"
TCP <- inegi_series(inegi_id, token, database = "BIE")
TCP <- ts_invert(TCP)
TCP <- ts(TCP, frequency = 12, start = c(2005,1))
ts.plot(TCP)

Verificar que sea estacionaria Contraste Aumentado de Dickey-Fuller

library(tseries)

## Warning: package 'tseries' was built under R version 4.3.3

adf.test(TCP, k=12)

## Warning in adf.test(TCP, k = 12): p-value smaller than printed p-value

## 
##  Augmented Dickey-Fuller Test
## 
## data:  TCP
## Dickey-Fuller = -4.1521, Lag order = 12, p-value = 0.01
## alternative hypothesis: stationary

par(mfrow=c(1,2))
Acf(TCP)
Pacf(TCP)

ARIMA(p,d,q)=ARIMa(1,0,0)

library(forecast)
m1 <- Arima(TCP, order = c(1,0,0))
m1

## Series: TCP 
## ARIMA(1,0,0) with non-zero mean 
## 
## Coefficients:
##          ar1     mean
##       0.6286  22.5630
## s.e.  0.0505   0.1118
## 
## sigma^2 = 0.4195:  log likelihood = -233.58
## AIC=473.16   AICc=473.26   BIC=483.58

ts.plot(TCP, m1$fitted, col=1:2)

summary(auto.arima(TCP, ic="bic"))

## Series: TCP 
## ARIMA(0,1,2)(1,0,1)[12] 
## 
## Coefficients:
##           ma1      ma2    sar1     sma1
##       -0.3124  -0.4478  0.8255  -0.6774
## s.e.   0.0640   0.0754  0.0873   0.1087
## 
## sigma^2 = 0.3874:  log likelihood = -223.01
## AIC=456.02   AICc=456.28   BIC=473.36
## 
## Training set error measures:
##                       ME      RMSE       MAE        MPE     MAPE      MASE
## Training set -0.01923807 0.6158412 0.4374254 -0.1560695 1.973706 0.6313186
##                    ACF1
## Training set 0.02425903

4.1.2 Modelos MA

Ocupación, empleo y remuneraciones > Tasas de ocupación, desocupación y subocupación (resultados mensuales de la ENOE, 15 años y más) > Urbana, agregado de 32 ciudades > Tasas complementarias > Tasa de ocupación en el sector informal 1 (TOSI1) > Total Clave:444803

4.1.3 Modelos ARMA

Indicadores económicos de coyuntura > Balanza comercial de mercancías de México > Series originales > Saldo > Total Clave:897

library(inegiR)
inegi_id <- "897"
AT <- inegi_series(inegi_id, token, database = "BIE")
AT <- ts_invert(AT)
AT <- ts(AT, frequency = 12, start = c(2005,1))
ts.plot(AT)

Verificar que sea estacionaria Contraste Aumentado de Dickey-Fuller

library(tseries)
adf.test(AT, k=4)

## Warning in adf.test(AT, k = 4): p-value smaller than printed p-value

## 
##  Augmented Dickey-Fuller Test
## 
## data:  AT
## Dickey-Fuller = -6.4478, Lag order = 4, p-value = 0.01
## alternative hypothesis: stationary

library(forecast)
par(mfrow=c(1,2))
Acf(diff(AT))
Pacf(diff(AT))

ARIMA(p,d,q)=ARIMa(2,1,2)

library(forecast)
m1 <- Arima(AT, order = c(2,1,2))
m1

## Series: AT 
## ARIMA(2,1,2) 
## 
## Coefficients:
##           ar1      ar2      ma1      ma2
##       -0.4406  -0.0130  -0.3134  -0.3352
## s.e.   0.5125   0.0861   0.5099   0.4022
## 
## sigma^2 = 2023007:  log likelihood = -3513.41
## AIC=7036.82   AICc=7036.97   BIC=7056.84

4.2 Fenómenos no estacionarios

4.2.1 Modelos integrados

4.2.2 Modelos ARIMA

Indicadores económicos de coyuntura > Índices de precios > Paridades de poder de compra para el Producto Interno Bruto de los países de la OCDE > Producto interno bruto > En moneda nacional a precios corrientes > Estados Unidos Clave:366709

inegi_id <- "366709"
PEA <- inegi_series(inegi_id, token, database = "BIE")
PEA <- ts_invert(PEA)
PEA <- ts(PEA, frequency = 1, start = c(2005,1))
ts.plot(PEA)

Verificar que sea estacionaria Contraste Aumentado de Dickey-Fuller

library(tseries)
adf.test(PEA, k=2)

## 
##  Augmented Dickey-Fuller Test
## 
## data:  PEA
## Dickey-Fuller = -0.90883, Lag order = 2, p-value = 0.934
## alternative hypothesis: stationary

plot(diff(PEA))

adf.test(diff(PEA), k=2)

## 
##  Augmented Dickey-Fuller Test
## 
## data:  diff(PEA)
## Dickey-Fuller = -2.5561, Lag order = 2, p-value = 0.3605
## alternative hypothesis: stationary

adf.test(diff(diff(PEA)), k=2)

## 
##  Augmented Dickey-Fuller Test
## 
## data:  diff(diff(PEA))
## Dickey-Fuller = -4.0223, Lag order = 2, p-value = 0.02248
## alternative hypothesis: stationary

par(mfrow=c(1,2))
Acf(PEA)
Pacf(PEA)

ARIMA(p,d,q)=ARIMa(2,1,2)

library(forecast)
m1 <- Arima(PEA, order = c(2,1,2))
m1

## Series: PEA 
## ARIMA(2,1,2) 
## 
## Coefficients:
##          ar1     ar2      ma1      ma2
##       0.5618  0.4350  -0.8380  -0.0308
## s.e.  0.8387  0.8391   0.7497   0.7194
## 
## sigma^2 = 245226518374:  log likelihood = -362.54
## AIC=735.08   AICc=738.24   BIC=741.18

ts.plot(PEA, m1$fitted, col=1:2)

4.3 Fenómenos con estacionalidad

4.3.1 Modelos estacionarios con estacionalidad

4.3.1.1 Modelos SARMA

Indicadores económicos de coyuntura > Encuesta mensual de la industria manufacturera (EMIM). Serie 2018 > Series Originales > Índices de remuneraciones medias reales por hora trabajada > Índice > Índice de remuneraciones medias reales por hora trabajada del personal dependiente de la razón social, considerando utilidades repartidas > 313 Fabricación de insumos textiles y acabado de textiles Clave:702286

inegi_id <- "702286"
AT <- inegi_series(inegi_id, token, database = "BIE")
AT <- ts_invert(AT)
AT <- ts(AT, frequency = 12, start = c(2005,1))
ts.plot(AT)

Verificar que sea estacionaria Contraste Aumentado de Dickey-Fuller

library(tseries)
adf.test(AT, k=12)

## 
##  Augmented Dickey-Fuller Test
## 
## data:  AT
## Dickey-Fuller = -3.4554, Lag order = 12, p-value = 0.04809
## alternative hypothesis: stationary

par(mfrow=c(1,2))
Acf(AT)
Pacf(AT)

4.3.2 Modelos no estacionarios con estacionalidad

4.3.2.1 Modelos SARIMA

Indicadores económicos de coyuntura > Producto interno bruto trimestral, base 2018 > Series Originales > Valores a precios de 2018 > Actividades primarias > 11 Agricultura, cría y explotación de animales, aprovechamiento forestal, pesca y caza Clave:735882 Paso 1 Verificar que sea estacionario en niveles

inegi_id <- "735882"
PEA <- inegi_series(inegi_id, token, database = "BIE")
PEA <- ts_invert(PEA)
PEA <- ts(PEA, frequency = 12, start = c(2005,1))
ts.plot(PEA)

Verificar que sea estacionaria Contraste Aumentado de Dickey-Fuller

library(tseries)
adf.test(PEA, k=12)

## 
##  Augmented Dickey-Fuller Test
## 
## data:  PEA
## Dickey-Fuller = -1.9733, Lag order = 12, p-value = 0.5874
## alternative hypothesis: stationary

plot(diff(PEA))

adf.test(diff(PEA), k=12)

## Warning in adf.test(diff(PEA), k = 12): p-value smaller than printed p-value

## 
##  Augmented Dickey-Fuller Test
## 
## data:  diff(PEA)
## Dickey-Fuller = -4.1943, Lag order = 12, p-value = 0.01
## alternative hypothesis: stationary

par(mfrow=c(1,2))
Acf(diff(PEA))
Pacf(diff(PEA))

ARIMA(p,d,q)=ARIMa(2,1,2)

library(forecast)
m1 <- Arima(PEA, order = c(2,1,2))
m1

## Series: PEA 
## ARIMA(2,1,2) 
## 
## Coefficients:
##          ar1     ar2      ma1     ma2
##       0.2335  0.7627  -1.7337  0.7396
## s.e.  0.0496  0.0500   0.0559  0.0557
## 
## sigma^2 = 1748135974:  log likelihood = -2146.49
## AIC=4302.97   AICc=4303.32   BIC=4318.88

ts.plot(PEA, m1$fitted, col=1:2)

summary(auto.arima(PEA, ic="bic"))

## Series: PEA 
## ARIMA(1,1,2)(0,1,1)[12] 
## 
## Coefficients:
##           ar1     ma1      ma2     sma1
##       -0.9746  0.0928  -0.5984  -0.6109
## s.e.   0.0321  0.0892   0.1029   0.0774
## 
## sigma^2 = 682286200:  log likelihood = -1924.99
## AIC=3859.98   AICc=3860.35   BIC=3875.54
## 
## Training set error measures:
##                   ME     RMSE      MAE         MPE    MAPE      MASE
## Training set 401.279 24849.31 19062.05 -0.07379344 3.19177 0.5409687
##                     ACF1
## Training set 0.003021428

inegi_id <- "736412"
ener <- inegi_series(inegi_id, token, database = "BIE")
tail(ener)

##           date date_shortcut   values notes
## 376 1993-06-01            M6 32.30856    NA
## 377 1993-05-01            M5 32.42110    NA
## 378 1993-04-01            M4 31.47853    NA
## 379 1993-03-01            M3 30.87515    NA
## 380 1993-02-01            M2 30.53460    NA
## 381 1993-01-01            M1 31.14023    NA

ener <- ts_invert(ener)
ener <- ts(ener, frequency = 12, start = c(1993,1)) 
ts.plot(ener)

Paso 1) análisis gráfico paso 2) PRUERBA de Raíces unitarias ADF y correlogramas Tiene tendencia, estacionalidad

adf.test(ener)

## Warning in adf.test(ener): p-value greater than printed p-value

## 
##  Augmented Dickey-Fuller Test
## 
## data:  ener
## Dickey-Fuller = 0.64895, Lag order = 7, p-value = 0.99
## alternative hypothesis: stationary

Por lo tanto, la variable ener no es estacionaria

Vamos a inspeccionar los correlogramas

Acf(ener)

Pacf(ener)

paso 3) aplicamos primeras diferencias

dener <- diff(ener)
ts.plot(dener)

adf.test(dener)

## Warning in adf.test(dener): p-value smaller than printed p-value

## 
##  Augmented Dickey-Fuller Test
## 
## data:  dener
## Dickey-Fuller = -11.314, Lag order = 7, p-value = 0.01
## alternative hypothesis: stationary

Acf(dener)

Pacf(dener)

Modelar parte no estacional

m1 <-Arima(ener, order = c(2,1,2))
summary(m1)

## Series: ener 
## ARIMA(2,1,2) 
## 
## Coefficients:
##          ar1     ar2     ma1     ma2
##       0.1065  0.1196  0.0408  0.3533
## s.e.  0.2190  0.1026  0.2240  0.0673
## 
## sigma^2 = 7.101:  log likelihood = -909.88
## AIC=1829.77   AICc=1829.93   BIC=1849.47
## 
## Training set error measures:
##                      ME    RMSE      MAE       MPE     MAPE      MASE
## Training set 0.06092265 2.64727 1.844564 0.1101307 2.631435 0.4375425
##                     ACF1
## Training set 0.006414525

ts.plot(ener, m1$fitted, col=1:2)

m2 <-Arima(ener, order = c(2,1,2), seasonal = c(1,1,0), lambda = 0)
summary(m2)

## Series: ener 
## ARIMA(2,1,2)(1,1,0)[12] 
## Box Cox transformation: lambda= 0 
## 
## Coefficients:
##           ar1      ar2     ma1     ma2     sar1
##       -0.2767  -0.0887  0.1794  0.2591  -0.4928
## s.e.   0.4498   0.3082  0.4408  0.2563   0.0456
## 
## sigma^2 = 0.0008792:  log likelihood = 773.45
## AIC=-1534.9   AICc=-1534.67   BIC=-1511.45
## 
## Training set error measures:
##                       ME     RMSE      MAE          MPE     MAPE      MASE
## Training set -0.02433139 2.040177 1.320867 -0.004510883 1.939994 0.3133182
##                     ACF1
## Training set -0.09566372

ts.plot(ener, m2$fitted, col=1:2)

plot(forecast(m2))

inegi_id <- "746097"
pib <- inegi_series(inegi_id, token, database = "BIE")
tail(pib)

##    date date_shortcut   values notes
## 15 2008          Year 20442062    NA
## 16 2007          Year 20251027    NA
## 17 2006          Year 19838804    NA
## 18 2005          Year 18929251    NA
## 19 2004          Year 18537508    NA
## 20 2003          Year 17899318    NA

pib <- ts_invert(pib)
pib <- ts(pib, frequency = 1, start = c(2003)) 
ts.plot(pib)

adf.test(pib)

## 
##  Augmented Dickey-Fuller Test
## 
## data:  pib
## Dickey-Fuller = -2.4387, Lag order = 2, p-value = 0.4053
## alternative hypothesis: stationary

dpib <- diff(pib)
ts.plot(dpib)

adf.test(dpib, k=2)

## 
##  Augmented Dickey-Fuller Test
## 
## data:  dpib
## Dickey-Fuller = -3.1207, Lag order = 2, p-value = 0.1454
## alternative hypothesis: stationary

ddpib<- diff(dpib)
ts.plot(ddpib)

adf.test(ddpib, k=2)

## 
##  Augmented Dickey-Fuller Test
## 
## data:  ddpib
## Dickey-Fuller = -3.7084, Lag order = 2, p-value = 0.04225
## alternative hypothesis: stationary

Acf(pib)

Pacf(pib)

par(mfrow=c(1,2))
Acf(pib)
Pacf(pib)

ARIMA(1,0,0) ARIMa(1,0,3) ARIMa(1,1,0) ARIMa(1,2,0)} ARIMA(1,1,3) ARIMA(1,2,3)

library(broom)
m1<- Arima(pib, order=c(1,0,0))
gm1 <- glance(m1)
AIC1 <- gm1$AIC
ts.plot(pib,m1$fitted, col=1:2)

m2<- Arima(pib, order=c(1,0,3))
gm2 <- glance(m2)
AIC2 <- gm2$AIC
m3<- Arima(pib, order=c(1,1,0))
gm3 <- glance(m3)
AIC3 <- gm3$AIC
m4<- Arima(pib, order=c(1,1,3))
gm4 <- glance(m4)
AIC4 <- gm4$AIC
m5<- Arima(pib, order=c(1,2,0))
gm5 <- glance(m5)
AIC5 <- gm5$AIC
m6<- Arima(pib, order=c(1,2,3))
gm6 <- glance(m6)
AIC6 <- gm6$AIC
AICS<-data.frame(c("m1", "m2", "m3", "m4", "m5", "m6"),
           c(AIC1, AIC2, AIC3, AIC4, AIC5, AIC6))
colnames(AICS) <- c("modelo", "AIC")
AICS

##   modelo      AIC
## 1     m1 609.9071
## 2     m2 615.7220
## 3     m3 575.7077
## 4     m4 579.8303
## 5     m5 553.2977
## 6     m6 550.3682

ts.plot(pib,m1$fitted,
        m2$fitted,
        m3$fitted,
        m4$fitted,
        m5$fitted,
        m6$fitted,
        col=1:7)

655561 Ocupación, empleo y remuneraciones > Indicadores de competitividad laboral > Salarios en México por subsector de actividad en la industria manufacturera > Mensual > 336 Fabricación de equipo de transporte

inegi_id <- "655561" 
trans <- inegi_series(inegi_id, token, database = "BIE")
tail(trans)

##          date date_shortcut values notes
## 76 2018-06-01            M6    2.6    NA
## 77 2018-05-01            M5    2.6    NA
## 78 2018-04-01            M4    2.6    NA
## 79 2018-03-01            M3    2.8    NA
## 80 2018-02-01            M2    2.7    NA
## 81 2018-01-01            M1    2.6    NA

trans <- ts_invert(trans)
trans <- ts(trans, frequency = 12, start = c(2018,1)) 
ts.plot(trans)

adf.test(trans, k=12)

## 
##  Augmented Dickey-Fuller Test
## 
## data:  trans
## Dickey-Fuller = -1.384, Lag order = 12, p-value = 0.828
## alternative hypothesis: stationary

dtrans<- diff(trans)
ts.plot(dtrans)

adf.test(dtrans, k=12)

## 
##  Augmented Dickey-Fuller Test
## 
## data:  dtrans
## Dickey-Fuller = -2.8475, Lag order = 12, p-value = 0.2289
## alternative hypothesis: stationary

adf.test(diff(dtrans),k=12)

## Warning in adf.test(diff(dtrans), k = 12): p-value smaller than printed p-value

## 
##  Augmented Dickey-Fuller Test
## 
## data:  diff(dtrans)
## Dickey-Fuller = -4.8252, Lag order = 12, p-value = 0.01
## alternative hypothesis: stationary

ddtrans<- diff(dtrans)
ts.plot(ddtrans)

par(mfrow=c(1,2))
Acf(ddtrans)
Pacf(ddtrans)

SARIMA(4,2,1)(1,0,2) SARIMA(0,2,1)(1,0,2) SARIMA(3,2,1)(1,0,2) SARIMA(2,2,1)(1,0,2)

m1 <- Arima(trans, order = c(0,2,1), seasonal = c(1,1,2))
m2 <- Arima(trans, order = c(3,2,1), seasonal = c(1,1,2))
m3 <- Arima(trans, order = c(2,2,1), seasonal = c(1,1,2))

ts.plot(trans, m1$fitted, m2$fitted,m3$fitted, col=1:4)

gm1 <- glance(m1)
AIC1 <- gm1$AIC
gm2 <- glance(m2)
AIC2 <- gm2$AIC
gm3 <- glance(m3)
AIC3 <- gm3$AIC
c(AIC1, AIC2, AIC3)

## [1] 135.3360 125.8999 125.1127

plot(forecast(m2))

5 Final Words

We have finished a nice book.

References

---

1. su ubicación específica es Indicadores económicos de coyuntura > Índices de precios > Índice nacional de precios al consumidor > Mensual > Índice↩︎